#@title
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
#@title
import numpy as np
import scipy
from sklearn.mixture import GaussianMixture
import matplotlib.pyplot as plt
import seaborn as sns
import pandas as pd
from scipy.stats import multivariate_normal
import matplotlib.colors as mcolors
from scipy import integrate
from scipy.optimize import minimize
import plotly.graph_objs as go
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
import plotly.figure_factory as ff
init_notebook_mode(connected=True)
from scipy.optimize import curve_fit
from plotly import tools
def configure_plotly_browser_state():
import IPython
display(IPython.core.display.HTML('''
<script src="/static/components/requirejs/require.js"></script>
<script>
requirejs.config({
paths: {
base: '/static/base',
plotly: 'https://cdn.plot.ly/plotly-latest.min.js?noext',
},
});
</script>
'''))
# np.inf distancia del maximo
def norma(x,orden=2):
return np.linalg.norm(x,axis=1,ord=orden)
# se guardan los registros de distancia y tiempo en un Dataframe de Pandas
registros=[[1.22, 2.62],[3.19, 4.26],[3.48,8.35],[4.11,7.00],[5.10,6.52],[3.09,5.49],[3.39,3.51],[4.96,7.64],[4.13,6.52],[1.64,3.09],[1.75,2.46],[3.23,3.89],[2.95,5.02],[3.07,5.49],[1.30,1.73],[4.26,6.82],[0.76,1.14],[4.40,5.53],[2.52,4.56],[2.42,4.30],[1.66,2.90],[2.96,3.55],[1.84,3.19],[0,0]]
registros=pd.DataFrame(registros,columns=['Distancia','Tiempo']).sort_values('Distancia').reset_index(drop=True)
#@title
# Para graficar los registros como nube de puntos
configure_plotly_browser_state()
#@title
trace = go.Scatter(
x = registros.Distancia[1:],
y = registros.Tiempo[1:],
mode = 'markers',
name = 'Registros'
)
point = go.Scatter(
x = [0],
y = [0],
mode = 'markers',
name = 'Origen'
)
layout = go.Layout(
title='Registros Distancia Vs. Tiempo',
font=dict(family='Courier New, monospace', size=18, color='#7f7f7f'),
xaxis=dict(
title='Distancia',
autorange=True,
titlefont=dict(
family='Courier New, monospace',
size=18,
color='#7f7f7f'
)
),
yaxis=dict(
title='Tiempo',
autorange=True,
titlefont=dict(
family='Courier New, monospace',
size=18,
color='#7f7f7f'
)
)
)
data = [trace,point]
fig = go.Figure(data=data, layout=layout)
# Plot and embed in ipython notebook!
iplot(fig)
# puntos para evaluarlos en las aproximaciones de la función t(x)
# se agragan al DataFrame ya que este almacena las evaluaciones y errores
xs = np.sort(list(np.arange(0, norma([[6,6]]), 0.1))+list(registros.Distancia))
for i in xs:
registros=registros.append({'Distancia':i},ignore_index=True)
registros=registros.drop_duplicates('Distancia').sort_values('Distancia').reset_index(drop=True)
originales=registros[['Distancia','Tiempo']].dropna()
# Calculo spline cubico
cubic_spline=scipy.interpolate.CubicSpline(originales.Distancia, originales.Tiempo, axis=0, bc_type='not-a-knot', extrapolate=None)
# Evaluación spline cubico en los puntos anteriormente dichos y calculo de error con los puntos originales
registros['Spline']=cubic_spline(registros.Distancia)
registros['Spline_error']=abs(registros.Tiempo-registros.Spline)
## calcula los coeficientes de menor a mayor grado y los valores evaluados de un ajuste
## por minimos cuadrados dado las listas x, y de igual tamaño y el grado del polinomio
## que se quiere ajustar
def polinomio_ajuste(x,y,grado):
columns=grado#+1
filas = len(x)
matriz=np.ones((filas,columns))
for i in range(columns):
matriz[:,i]=np.power(x,i+1)
#matriz[:,i]=np.power(x,i) para considerar coeficiente independiente
#matriz[:,i]=np.power(x,i+1) para no considerar coeficiente independiente
y=np.array(y).reshape((len(y),1))
coef=np.dot(np.dot(np.linalg.inv(np.dot(matriz.T,matriz)),matriz.T),y)
y_approx=np.dot(matriz,coef)
return coef, y_approx
## Evalua una función polinomica para una lista de valores, el número de coeficientes
## determina el grado del polinomio
def evaluar_polinomio(x,coeficientes):
return np.dot(np.array([np.power(x,i+1) for i in range(len(coeficientes))]).T,coeficientes).flatten()
#np.power(x,i) para considerar coeficiente independiente
#np.power(x,i+1) para no considerar coeficiente independiente
# Se calculan los coeficientes y evaluaciones para los polinomios de grado 1,2 y 3
## y se almacenan en la tabla de registro junto con su errores en valor absoluto
polinomios=[]
for i,j in zip([1,2,3],['Lineal','Cuadratico','Cúbico']):
coef,t_approx=polinomio_ajuste(originales.Distancia,originales.Tiempo,i)
polinomios.append(coef)
registros[j]=evaluar_polinomio(registros.Distancia,coef)
registros[j+'_error']=abs(registros[j]-registros.Tiempo)
#@title
# grafica aproximaciones
configure_plotly_browser_state()
aproximaciones=[]
for i in ['Spline','Lineal','Cuadratico','Cúbico']:
grafica=go.Scatter(
x = registros.Distancia,
y = registros[i],
mode = 'lines',
name = i,
error_y=dict(
type='data',
array=registros[i+'_error'],
visible=True
)
)
aproximaciones.append(grafica)
data = data+aproximaciones
fig = go.Figure(data=data, layout=layout)
# Plot and embed in ipython notebook!
iplot(fig)
#@title
# Grafica errores de aproximación en los registros originales
configure_plotly_browser_state()
Errores=pd.DataFrame(registros[['Spline_error','Lineal_error','Cuadratico_error','Cúbico_error']].sum(),columns=['Error Estimación'])
trace = go.Bar(
x=Errores['Error Estimación']+0.01,
y=list(Errores.index),
orientation='h',
text=list(Errores.index),
textposition = 'auto',
marker=dict(
color='rgb(158,202,225)',
line=dict(
color='rgb(8,48,107)',
width=1),
),
)
layout = go.Layout(
title='Errores de Estimación',
font=dict(family='Courier New, monospace', size=18, color='#7f7f7f'),
xaxis=dict(
title='Error',
autorange=True,
titlefont=dict(
family='Courier New, monospace',
size=18,
color='#7f7f7f'
)
),
yaxis=dict(
autorange=True,
showgrid=False,
zeroline=False,
showline=False,
ticks='',
showticklabels=False
)
)
data = [trace]
fig = go.Figure(data=data, layout=layout)
# Plot and embed in ipython notebook!
iplot(fig)
## inicializar variable que almacena las soluciones y escogencia de parametros del modelo
soluciones={}
grado_funcion=2#Errores.reset_index(drop=True)[1:].idxmin().values[0]
distance_order=2 #np.inf
# funcion del tiempo necesario desde el hospital a una emergencia
def t(emergencia,hospital,coeficientes=polinomios[grado_funcion-1],distance_order=2):
return evaluar_polinomio(norma(np.array(emergencia)-np.array(hospital),distance_order),coeficientes)
## frecuencia llamadas en las casillas de la ciudad
frequency = np.array([[ 0,2,3,1,1,1],
[10,6,3,1,3,1],
[ 8,5,2,1,0,0],
[ 5,3,3,0,1,2],
[ 2,1,1,2,3,2],
[ 3,0,1,4,2,1]])
## calculo puntos centrales de las casillas
means=[]
for row in range(6):
for col in range(6):
means.append([col+0.5, row+0.5])
means=np.array(means)
X,Y=np.meshgrid(np.arange(7),np.arange(7))
fig, (ax1, ax2) = plt.subplots(1, 2,figsize=(20,8))
ax1.imshow(plt.imread("img.png"))
ax1.axis('off')
im=ax2.pcolormesh(X,Y,frequency,cmap='viridis',edgecolors='face',alpha=0.5,snap=True)
ax2.scatter(means[:,0],means[:,1],c='r')
fig.colorbar(im, ax=ax2)
plt.show()
# creación malla y puntos donde se evalua la función objetivo
step=0.025
X, Y = np.meshgrid(np.arange(0, 6., step), np.arange(0, 6., step))
points=[i for i in zip(X.flatten(),Y.flatten())]
## restricciones de la optimización, el minimo debe quedar dentro de la ciudad
def c1(x):
return x[0]
def c2(x):
return x[1]
def c3(x):
return 6- x[0]
def c4(x):
return 6- x[1]
cons = [{'type':'ineq', 'fun': c1},
{'type':'ineq', 'fun': c2},
{'type':'ineq', 'fun': c3},
{'type':'ineq', 'fun': c4}]
## metodos de optimizacion
metodos=['Nelder-Mead','Powell','COBYLA']
# definicion función objetivo modelo discreto
def discreto(hospital,frequency=frequency,coeficientes=polinomios[grado_funcion-1],distance_order=2):
return sum(t(means,hospital,coeficientes,distance_order)*frequency.flatten())
## tiempos de calculo de la solución para los metodos basados en simplex
time_metodos=[]
for i in metodos:
a=%timeit -o -q minimize(discreto, [2,2], method = i,args=(frequency,polinomios[grado_funcion-1],distance_order),tol=1e-10)
time_metodos.append([a.best])
time_metodos=np.array(time_metodos)
#@title
# grafica tiempos de calculo
configure_plotly_browser_state()
data = [
go.Bar( x=time_metodos[:,0], y=metodos,
orientation='h' ,
text=metodos,
textposition = 'auto',
marker=dict(
color='rgb(158,202,225)',
line=dict(
color='rgb(8,48,107)',
width=1),
),
)
]
layout = go.Layout(
title='Tiempos Convergencia Métodos Optimización',
font=dict(family='Courier New, monospace', size=18, color='#7f7f7f'),
xaxis=dict(
title='Tiempo (segundos)',
autorange=True,
titlefont=dict(
family='Courier New, monospace',
size=18,
color='#7f7f7f'
)
),
yaxis=dict(
autorange=True,
showgrid=False,
zeroline=False,
showline=False,
ticks='',
showticklabels=False
)
)
fig = go.Figure(data=data, layout=layout)
iplot(fig)
# optimizacion con algoritmo COBYLA
sol_d=minimize(discreto, [2,2], method = 'COBYLA',args=(frequency,polinomios[grado_funcion-1],distance_order),tol=1e-10)
soluciones['Discreto']=sol_d.x
print("Salida Optimización")
print(sol_d)
#@title
# grafica función objetivo
configure_plotly_browser_state()
h_disc=[discreto(i,frequency=frequency,coeficientes=polinomios[grado_funcion-1],distance_order=2) for i in points]
fig = tools.make_subplots(rows=1, cols=2,
specs=[[{'is_3d': False}, {'is_3d': True}]],print_grid=False)
# adding surfaces to subplots.
fig.append_trace(dict(type='contour', x=X[0], y=X[0], z=np.array(h_disc).reshape(X.shape), colorscale='Viridis',
showscale=True), 1, 1)
fig.append_trace(dict(type='scatter', x=[sol_d.x[0]], y=[sol_d.x[1]],mode = 'markers'), 1, 1)
fig.append_trace(dict(type='surface', x=X[0], y=X[0], z=np.array(h_disc).reshape(X.shape), colorscale='Viridis',
showscale=False), 1, 2)
fig['layout'].update(title='Evaluación H(x,y)',font=dict(family='Courier New, monospace', size=18, color='#7f7f7f'))
iplot(fig)
Para un punto $p_c$ con componentes $(x_c, y_c)$ en la $k$-cuadricula la función de densidad de probabilidad (Relacionada entonces con la frecuencia) es:
\begin{equation} \phi_c(p|\mu_k, \Sigma_k) = \frac{exp(-\frac{1}{2}((p-\mu_k)^T\Sigma_k^{-1}(p-\mu_k)))}{\sqrt{(2\pi)^2 |\Sigma_k|}} \end{equation}donde
$${\mu }={\begin{pmatrix}\mu _{x}\\\mu _{y}\end{pmatrix}},\quad {\boldsymbol \Sigma }={\begin{pmatrix}\sigma _{x}^{2}& 0\\0&\sigma _{y}^{2}\end{pmatrix}}. $$En nuestro modelo, asumimos que el $95\%$ de los datos se encuentran en la casilla, es decir $2\sigma = 0.5km$.
Para extrapolar el método de esta casilla a las 36 casillas en el mapa, usamos las mezclas gaussianas. Donde la función de densidad de probabilidad para cada punto de la región $p$ con coordenadas $(x,y)$ esta dado por:
\begin{equation} \phi(x) = \sum_{k=1}^{36} \pi_k \phi_k(p | \mu_k, \Sigma_k) \end{equation}Los parametros $\pi_k$ son llamados coeficientes de mixtura, donde se sabe que:
$$ \sum_{k=1}^{36} \pi_k =1 ,\quad 0 \leq \pi_k \leq 1. $$Si normalizamos los valores de las frecuencias de llamadas por cuadricula, obtenemos los valores de $\pi_k$, es decir.
$$ \pi_k = \frac{f(x_i, y_y) }{ \sum_{k=1}^{36} f(x_i, y_y) } $$Donde la nueva función a optimizar seria:
$$\underset{(x_0, y_0) \in A}{\text{minimizar}} \iint{\phi(x,y) t\left( d((x_0,y_0),(x,y)) \right) dA}$$# normalizacion frecuencias y funcion phi
pi=frequency.flatten()/frequency.sum()
def fdp(x):
return 8/np.pi*sum(pi*np.exp(-8*np.sum(np.power(means-x,2),axis=1)))
step=0.025
X, Y = np.meshgrid(np.arange(0, 6., step), np.arange(0, 6., step))
points=[i for i in zip(X.flatten(),Y.flatten())]
# evaluación función phi en los puntos de la malla
fdp_values=np.array([fdp(i) for i in points])
# Centro de masa
cm=sum([i*j for i,j in zip(fdp_values,np.array(points))])/sum(fdp_values)
soluciones['Centro Masa']=cm
#@title
#grafica funcion phi
configure_plotly_browser_state()
Z=fdp_values.reshape(X.shape)
fig = tools.make_subplots(rows=1, cols=2,
specs=[[{'is_3d': False}, {'is_3d': True}]],print_grid=False)
# adding surfaces to subplots.
fig.append_trace(dict(type='contour', x=X[0], y=X[0], z=Z, colorscale='Viridis',
showscale=True), 1, 1)
fig.append_trace(dict(type='surface', x=X[0], y=X[0], z=Z, colorscale='Viridis',
showscale=False), 1, 2)
fig['layout'].update(title='Densidad Probabilidad Llamadas y Curvas de nivel',font=dict(family='Courier New, monospace', size=18, color='#7f7f7f'))
iplot(fig)
## devuelve la funcion del integrando que depede de la ubicación del hospital
def parameters(hospital=[3,3],coeficientes=polinomios[grado_funcion-1],distance_order=distance_order):
# funcion integrando en la definicion del modelo continuo
def integrando(y,x,hospital=hospital,coeficientes=coeficientes,distance_order=distance_order):
return (fdp([x,y])*t([[x,y]],hospital,coeficientes=coeficientes,distance_order=distance_order))[0]
return integrando
## calculo integral dada una ubicación del hospital
def integrar(x,coeficientes=polinomios[grado_funcion-1],distance_order=2):
return integrate.dblquad(parameters(hospital=x,coeficientes=coeficientes,distance_order=distance_order), 0, 6, 0, 6)[0]
# optimizacion función objetivo modelo continuo
sol=minimize(integrar, [2,2],method = 'COBYLA',args=(polinomios[grado_funcion-1],distance_order),constraints=cons,tol=1e-10)
soluciones['Continuo']=sol.x
sol
"""
step=0.1
X, Y = np.meshgrid(np.arange(0, 6., step), np.arange(0, 6., step))
points=[i for i in zip(X.flatten(),Y.flatten())]
h_con=[integrar(i,coeficientes=polinomios[grado_funcion-1],distance_order=distance_order) for i in points]
fig = tools.make_subplots(rows=1, cols=2,
specs=[[{'is_3d': False}, {'is_3d': True}]],print_grid=False)
# adding surfaces to subplots.
fig.append_trace(dict(type='contour', x=X[0], y=X[0], z=np.array(h_con).reshape(X.shape), colorscale='Viridis',
showscale=True), 1, 1)
fig.append_trace(dict(type='scatter', x=[sol.x[0]], y=[sol.x[1]],mode = 'markers'), 1, 1)
fig.append_trace(dict(type='surface', x=X[0], y=X[0], z=np.array(h_con).reshape(X.shape), colorscale='Viridis',
showscale=False), 1, 2)
fig['layout'].update(title='Evaluación H(x,y) continuo')
iplot(fig)
"""
#@title
configure_plotly_browser_state()
Z=fdp_values.reshape(X.shape)
trace0 = go.Contour(
x=X[0],
y=X[0],
z=Z,
colorscale='Viridis',
showscale=True
)
trace1 = go.Scatter(
x=np.array(list(soluciones.values()))[:,0],
y=np.array(list(soluciones.values()))[:,1],
mode="markers+text",
text=list(soluciones.keys()),
textposition=["bottom left","top center","middle right"],
textfont=dict(
family="sans serif",
size=25,
color="red"
),
marker = dict(
size = 10,
color = np.random.randn(500),
#colorscale='Viridis',
#showscale=True,
line = dict(
width = 2,
)
))
data = [trace0,trace1]
layout = go.Layout(
title='Soluciones Encontradas',
font=dict(family='Courier New, monospace', size=18),
xaxis=dict(
title='X',
autorange=True,
titlefont=dict(
family='Courier New, monospace',
size=18
)
),
yaxis=dict(
title='Y',
autorange=True,
titlefont=dict(
family='Courier New, monospace',
size=18
)
)
)
# Plot and embed in ipython notebook!
fig = go.Figure(data=data, layout=layout)
iplot(fig)
step=0.025
X, Y = np.meshgrid(np.arange(0, 6., step), np.arange(0, 6., step))
points=[i for i in zip(X.flatten(),Y.flatten())]
# Rendimiento modelos
n_puntos_prueba=10000
puntos_prueba_uniforme=[points[np.random.choice(len(points))] for i in range(n_puntos_prueba)]
index=np.random.choice(np.arange(len(points)), n_puntos_prueba, p=fdp_values/sum(fdp_values))
puntos_prueba_fdp=[points[i] for i in index]
uniforme=np.array([])
means_uniforme={}
for h in soluciones:
arreglo=np.array(t(puntos_prueba_uniforme,soluciones[h],coeficientes=polinomios[grado_funcion-1],distance_order=distance_order))
means_uniforme[h]=np.mean(arreglo)
uniforme=np.concatenate((uniforme,arreglo),axis=0)
fdp_sample=np.array([])
means_fdp={}
for h in soluciones:
arreglo=np.array(t(puntos_prueba_fdp,soluciones[h],coeficientes=polinomios[grado_funcion-1],distance_order=distance_order))
means_fdp[h]=np.mean(arreglo)
fdp_sample=np.concatenate((fdp_sample,arreglo),axis=0)
annotations=[
dict(
x=-0.2+j*1,
y=means_uniforme[i],
xref='x',
yref='y',
text=format(means_uniforme[i], '.4f'),
showarrow=True,
#arrowhead=7,
ax=-60,
ay=0
) for i,j in zip(means_uniforme,np.arange(3))
]
annotations+=[
dict(
x=0.2+j*1,
y=means_fdp[i],
xref='x',
yref='y',
text=format(means_fdp[i], '.4f'),
showarrow=True,
#arrowhead=7,
ax=60,
ay=0
) for i,j in zip(means_fdp,np.arange(3))
]
#@title
configure_plotly_browser_state()
fig = {
"data": [
{
"type": 'violin',
"x": np.array([[i]*10000 for i in list(soluciones.keys())]).flatten(),
"y": uniforme,
"legendgroup": 'Muestreo Uniforme',
"scalegroup": 'Muestreo Uniforme',
"name": 'Muestreo Uniforme',
"box": {
"visible": False
},
"meanline": {
"visible": True
},
"line": {
"color": 'blue'
}
},
{
"type": 'violin',
"x": np.array([[i]*10000 for i in list(soluciones.keys())]).flatten(),
"y": fdp_sample,
"legendgroup": 'Muestreo fdp',
"scalegroup": 'Muestreo fdp',
"name": 'Muestreo fdp',
"box": {
"visible": False
},
"meanline": {
"visible": True
},
"line": {
"color": 'red'
}
}
],
"layout" : {
"yaxis": {
"autorange":True,
"zeroline": False,
},
"violinmode": "group",
"legend":dict(orientation="h"),
"title":'Evaluación rendimiento modelos',
"font":dict(family='Courier New, monospace', size=18, color='#7f7f7f'),
"annotations":annotations
}
}
iplot(fig)
configure_plotly_browser_state()
orders=[1,2,7,np.inf]
titles=[]
for i in orders:
for j in range(len(polinomios)):
titles.append("l"+str(i)+" grado "+ str(j+1))
fig = tools.make_subplots(rows=4, cols=3,print_grid=False,subplot_titles=titles)
for i in range(len(orders)):
for j in range(len(polinomios)):
s=minimize(discreto, [2,2], method = 'COBYLA',args=(frequency,polinomios[j],orders[i]),tol=1e-10)
h_disc=[discreto(k,frequency=frequency,coeficientes=polinomios[j],distance_order=orders[i]) for k in points]
fig.append_trace(dict(type='contour', x=X[0], y=X[0], z=np.array(h_disc).reshape(X.shape), colorscale='Viridis',
showscale=False), i+1, j+1)
fig.append_trace(dict(type='scatter', x=[s.x[0]], y=[s.x[1]],mode = 'markers'), i+1, j+1)
fig['layout'].update(showlegend=False,height=900, width=1015,title='Experimentos',font=dict(family='Courier New, monospace', size=18, color='#7f7f7f'))
iplot(fig)